#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>


#define ITMAX 100 /*Maximum allowed number of iterations.*/
#define EPS 3.0e-7 /*Relative accuracy.*/
#define FPMIN 1.0e-30 /*Number near the smallest representable*/

static void gser(float *gamser, float a, float x, float *gln);
static void gcf(float *gammcf, float a, float x, float *gln);
static float gammq(float a, float x);
static double static_normalp(const double x);
static float gammp(float a, float x);

/*funcion de probabilidad de masa de la geometrica de parametro 0.4*/
double p_geom(double x){
	assert(x>=1.0);

	return 0.4*pow(0.6,x-1.0);
}

/*Probabilidad chi-cuadrado. Calcula P(T<=t)*/
double chi_squareq(double a, double x){
        return gammq(a/2.0, x/2.0);
}

/*Normal de parametro mu y sigma acumulada*/
double normalp(const double x, double mu, double sigma){
	return static_normalp((x - mu)/sigma);
}

/*Lognormal de parametro mu y sigma*/
double log_normalp(const double x, double mu, double sigma){
	return normalp(log(x), mu, sigma);
}


/*Gamma de parametro alpha y beta*/
double gammap(double x, double alpha, double beta){
	
	return gammp(alpha, (x/beta));
}

float gammln(float xx)
/*Returns the value ln[􀀀(xx)] for xx > 0.*/
{
        /*Internal arithmetic will be done in double precision, a nicety that you can omit if ve-gure
        accuracy is good enough.*/
        double x,y,tmp,ser;
        static double cof[6]={76.18009172947146,-86.50532032941677,
        24.01409824083091,-1.231739572450155,
        0.1208650973866179e-2,-0.5395239384953e-5};
        int j;
        y=x=xx;
        tmp=x+5.5;
        tmp -= (x+0.5)*log(tmp);
        ser=1.000000000190015;
        for (j=0;j<=5;j++) ser += cof[j]/++y;
                return -tmp+log(2.5066282746310005*ser/x);
}

/*Funcion imagen*/
double U1(double x){
	return x;
}

/*Funcion que ordena la muestra*/
void bizzard_sort(double *ma, unsigned int n){
	unsigned int i = 0, j = 0;
	
	for(j = 0 ; j < n ; j++){
		for(i = 1; i < n; i++){
			if(ma[i] < ma[i - 1]){
				double tmp = ma[i - 1];
				ma[i - 1] = ma[i];
				ma[i] = tmp;
			}
		}
	}
}


/*FUNCIONES AUXILIARES PARA CHI-CUADRADO y GAMMA****/
static float gammq(float a, float x)
/*Returns the incomplete gamma function Q(a; x) = 1 − P(a; x).*/
{
        void gcf(float *gammcf, float a, float x, float *gln);
        void gser(float *gamser, float a, float x, float *gln);
        void nrerror(char error_text[]);
        float gamser,gammcf,gln;
        if (x < 0.0 || a <= 0.0) printf("Invalid arguments in routine gammq\n");
        if (x < (a+1.0)) { /*Use the series representation*/
                gser(&gamser,a,x,&gln);
                return 1.0-gamser; /*and take its complement.*/
        } else { /*Use the continued fraction representation.*/
                gcf(&gammcf,a,x,&gln);
                return gammcf;
        }
}

static void gser(float *gamser, float a, float x, float *gln)
/*Returns the incomplete gamma function P(a; x) evaluated by its series representation as gamser.
Also returns ln 􀀀(a) as gln.*/
{
        float gammln(float xx);
        void nrerror(char error_text[]);
        int n;
        float sum,del,ap;
        *gln=gammln(a);
        if (x <= 0.0) {
        if (x < 0.0) printf("x less than 0 in routine gser\n");
                *gamser=0.0;
                return;
        } else {
                ap=a;
                del=sum=1.0/a;
                for (n=1;n<=ITMAX;n++) {
                        ++ap;
                        del *= x/ap;
                        sum += del;
                        if (fabs(del) < fabs(sum)*EPS) {
                                *gamser=sum*exp(-x+a*log(x)-(*gln));
                                return;
                        }
                }
                printf("a too large, ITMAX too small in routine gser\n");
                return;
        }
}

/*floating-point number.*/
static void gcf(float *gammcf, float a, float x, float *gln)
/*Returns the incomplete gamma function Q(a; x) evaluated by its continued fraction representation
as gammcf. Also returns ln 􀀀(a) as gln.*/
{
        float gammln(float xx);
        void nrerror(char error_text[]);
        int i;
        float an,b,c,d,del,h;
        *gln=gammln(a);
        b=x+1.0-a; /*Set up for evaluating continued fraction
        by modied Lentz's method (x5.2)
        with b0 = 0.*/
        c=1.0/FPMIN;
        d=1.0/b;
        h=d;
        for (i=1;i<=ITMAX;i++) { /*Iterate to convergence.*/
                an = -i*(i-a);
                b += 2.0;
                d=an*d+b;
                if (fabs(d) < FPMIN) d=FPMIN;
                c=b+an/c;
                if (fabs(c) < FPMIN) c=FPMIN;
                d=1.0/d;
                del=d*c;
                h *= del;
                if (fabs(del-1.0) < EPS) break;
        }
        if (i > ITMAX) printf("a too large, ITMAX too small in gcf\n");
        *gammcf=exp(-x+a*log(x)-(*gln))*h;/* Put factors in front.*/
}

/** **********/


/**FUNCION AUXILIAR PARA LA NORMAL*/
static double static_normalp(const double x){
	unsigned int i;
	double result = 0.0;
	double xsq;
	double xnum;
	double xden;
	
	const double a[5] = {
		2.2352520354606839287,
	        161.02823106855587881,
	        1067.6894854603709582,
	        18154.981253343561249,
	        0.065682337918207449113
	};
	const double b[4] = {
		47.20258190468824187,
	        976.09855173777669322,
	        10260.932208618978205,
	        45507.789335026729956
	};
	
	xsq = x * x;
	xnum = a[4] * xsq;
	xden = xsq;
	
	for (i = 0; i < 3; i++)
	{
		xnum = (xnum + a[i]) * xsq;
		xden = (xden + b[i]) * xsq;
	}
	
	result = x * (xnum + a[3]) / (xden + b[3]);
	
	
	return 0.5 + result;
}
/** *******/

/**FUNCION AUXILIAR PARA LA GAMMA*/
static float gammp(float a, float x)
/*Returns the incomplete gamma function P(a; x).*/
{
	void gcf(float *gammcf, float a, float x, float *gln);
	void gser(float *gamser, float a, float x, float *gln);
	float gamser,gammcf,gln;
	
	if (x < 0.0 || a <= 0.0) printf("Invalid arguments in routine gammp\n");
	
	if (x < (a + 1.0)) {/* Use the series representation.*/
		gser(&gamser,a,x,&gln);
		return gamser;
	} else { /*Use the continued fraction representation*/
		gcf(&gammcf, a, x, &gln);
		return 1.0 - gammcf; /*and take its complement.*/
	}
}
